Method and system for near-infrared fluorescence contrast-enhanced imaging with area illumination and area detection

ABSTRACT

According to one embodiment of the invention, a method for biomedical imaging includes directing time-varying excitation light at a surface area of a light scattering material, the material comprising a fluorescent target. Time-varying emission light from the fluorescent target is detected, substantially at a two-dimensional sensor surface, in response to the time-varying excitation light stimulating the fluorescent target. The time-varying emission light is filtered to reject excitation light re-emitted from the material. A three-dimensional image of the fluorescent target is generated based on the detection substantially at the sensor surface.

BACKGROUND OF THE INVENTION

Fluorescent contrast agents have been developed and/or employed thatenhance the optical detection of diseased tissues. These agents exciteand re-emit at near-infrared (NIR) wavelengths, which deeply penetratetissues. Localization of a contrast-enhanced target in three dimensionsnecessitates: (1) rapid and precise acquisition of fluorescencemeasurements; (2) accurate prediction of these measurements using anappropriate theoretical model of light propagation through tissues; and(3) tomographic reconstruction of model parameters, which includeoptical and fluorescent properties of tissues, using an optimizationroutine that minimizes the differences between the experimental andpredicted fluorescence data. Steady-state and time-resolved experimentalmeasurements of diffuse fluorescence have been used to localize acontrast-enhanced target. The steady-state measurements provideinformation regarding the spatial origin of fluorescent photons, and thetime-resolved measurements provide additional information regardingfluorescence decay kinetics, which may render functional informationabout the environment in which the fluorophores reside.

Experimental measurements of diffuse fluorescence acquired using pointillumination and point detection techniques have been incorporated intotomographic reconstruction algorithms. In order to adequately probe atissue volume using these techniques, illumination and detection mustoccur via fiber optic at several points about the tissue volume.Disadvantages with point illumination and point detection include (1)possible failure to sufficiently excite the contrast-enhanced target and(2) the production of sparse data sets, which make the underdeterminedinverse problem difficult to solve.

SUMMARY OF THE INVENTION

A method and system for near-infrared fluorescence contrasted-enhancedimaging with area illumination and area detection are provided. In oneembodiment, a method for biomedical imaging includes directingtime-varying excitation light at a surface area of a light scatteringmaterial, the material comprising a fluorescent target. Time-varyingemission light from the fluorescent target is detected, substantially ata two-dimensional sensor surface, in response to the time-varyingexcitation light stimulating the fluorescent target. The time-varyingemission light is filtered to reject excitation light re-emitted fromthe material. A three-dimensional image of the fluorescent target isgenerated based on the detection substantially at the sensor surface.

Embodiments of the invention provide a number of technical advantages.Embodiments of the invention may include all, some, or none of theseadvantages. An approach employed herein involves area illumination andarea detection on the same tissue-like surface. This approach overcomesthe disadvantages associated with point illumination and point detectiongiven that (1) fluorophores within a large tissue volume are moreefficiently excited via area illumination and (2) significantly moredata is provided via area detection. In addition, non-contact imaging ispossible. Some applications may arise in medical imaging of diseasedtissues and intraoperative 3D identification of diseased (and tumor)lesion margin for effective and complete resection of diseased tissue,among others.

Other advantages may be readily ascertainable by those skilled in theart from the following figures, description, and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the invention, and for furtherfeatures and advantages, reference is now made to the followingdescription, taken in conjunction with the accompanying drawings, inwhich:

FIG. 1 is a schematic illustrating a system for near-infrared (“NIR”)frequency-domain photon migration (“FDPM”) measurements according to oneembodiment of the present invention;

FIG. 2 is a schematic illustrating a system for detecting re-emittedexcitation light in accordance with one embodiment of the presentinvention;

FIG. 3 is a flow chart illustrating a method for tomographic fluorescentimaging in accordance with one embodiment of the present invention;

FIG. 4 illustrates one embodiment of generation step of FIG. 3; and

FIGS. 5A-C illustrate a reconstructed three-dimensional image of afluorescent target.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS OF THE INVENTION

FIG. 1 is a schematic illustrating a system for excitation areaillumination and area detection using near-infrared (“NIR”)frequency-domain photon migration (“FDPM”) imaging techniques accordingto one embodiment of the present invention. In particular, system 100generates three-dimensional (“3D”) images of fluorescent target 102within a light scattering material (e.g., tissue) based on emissionlight detected at a two-dimensional (“2D”) sensor surface. At a highlevel, system 100 delivers an intensity-modulated excitation light 104into the tissue of a body 106 to cause stimulated emission offluorescent target 102. Emission light 108 from the fluorescent target102 is detected by intensified charge-coupled camera (ICCD) 112 via gainmodulated image intensifier 114. A 3D image of fluorescent target 102 isgenerated based, at least in part, on excitation light 104 and emissionlight 108. System 100 provides several advantages over the prior artsuch as, for example, the ability to reconstruct a 3D image from datacollected at a 2D surface with contact. The present inventioncontemplates more, less, or different components than those shown insystem 100 of FIG. 1.

Referring to FIG. 1, system 100 includes a laser diode 116 coupled to alaser DC bias 118 and a reference frequency generator 120. As a result,laser diode 116 generates light source 122, intensity-modulatedexcitation light. In one embodiment, light source 122 comprises785-nanometer (“nm”) excitation light provided by a 70-milliWatt (“mW”)laser diode 116. Additionally, a temperature controller (notillustrated) may be include in system 100 to assist in maintaining theoptical power and the wavelength of laser diode 116; however, othersuitable wavelengths and sources may be used. Light source 122 maycomprise other suitable NIR light such as, for example, 700-nm to 900-nmlight. Light source 122 passes through lens 124A to produce an expandedbeam of excitation light 104 for illuminating a surface area 126 of body106. Surface area 126, in one embodiment, is approximately 9 centimeterssquared (cm²) in size. Surface area 126 may comprise other suitableareas ranging from 1 cm² to 100 cm². Surface area 126 is illuminatedwith excitation light 104 for stimulating fluorescent target 102.Fluorescent target 102 comprises a fluorescent contrast agent (e.g.,Indocyanine Green (“ICG”)) that emits light in response to excitationlight 104. In one embodiment, ICG is excited at 780 nm and emits at 830nm, has an extinction coefficient of 130,000 M⁻¹cm⁻¹, a fluorescentlifetime of 0.56 ns and a quantum efficiency of 0.016 for the 780/830 nmexcitation/emission wavelengths. After illuminating surface area 126,light is re-emitted from body 106.

Re-emitted light 110 includes emission light 108 and reflectedexcitation light. To minimize the excitation light delivered to imageintensifier 114, a band-pass filter 128 centered around the emissionwavelength and a band-rejection filter 130 centered around theexcitation wavelength are used in combination to substantially rejectthe excitation wavelength from re-emitted light 110 to provide emissionlight 108. Emission light 108 substantially comprises light emitted byfluorescent target 102 as a result of excitation light 104 stimulatingfluorescent target 102. Band-pass filter 128 and band-rejection filter130 may be coupled via lens 124B (e.g., 50-mm lens). In one embodiment,band-pass filter 128 and band-rejection filter 130 comprise an 830-nmimage quality bandpass interference filter and a 785-nm holographicband-rejection filter, respectively. Alternatively, system 100 may useany suitable combination of filters such as, for example band-passfilters, long-pass filters, holographic notch filters, band-rejectionfilters or any combination thereof, or any optical scheme to reject theexcitation wavelength enabling registration of the emission wavelength.After the combined filters substantially remove the excitationwavelength, emission light 108 is amplified by image intensifier 114.

In one embodiment, image intensifier 114 includes a photocathode facethat converts photons to electrons, a multi-channel plate (“MCP”) thatmultiplies the electronic signal by avalanche multiplication, and aphosphorescent screen that converts electrons into an optical image. Inone embodiment, image intensifier 114 is a fast intensifier of thevariety manufactured by Litton Electronics, Inc., which enablesmodulation by applying a DC bias via image intensifier bias 132 and anRF signal via an RF amplifier 134 between the photocathode and the MCP;however, other suitable image intensifiers may be utilized. In theillustrated embodiment, the intensifier gain of image intensifier 114and the intensity of excitation light 104 are phase-locked by a 10 MHzreference signal provided by oscillator 120. By modulating laser diode116 and image intensifier 114 at the same frequency, homodyning resultsin a steady-state image on the phosphor screen. The image from thephosphor screen is focused through a lens 124C (e.g., 105-mm lens) ontoa two-dimensional sensor surface 138 of camera 112, in order to providearea detection. As used herein, the term “two-dimensional sensorsurface” means a substantially two-dimensional surface determined by thesensors of camera 112. Sensor surface 138 may be flat, curved, bent, orarranged in any suitable surface. In one embodiment, sensor surface 138comprises a 16 cm² flat surface. Camera 112, in the illustratedembodiment, has a 512×512 array of CCD detectors configured to provide acorresponding pixelated image.

In operation, following each acquired image, a phase delay between imageintensifier 114 and laser diode 116 is induced by stepping the phase ofimage intensifier 114 between 0 and 360° relative to the phase of laserdiode 116 modulation. Preferably, control between oscillator 140 and theprocessor is obtained by a conventional general purpose interface bus(“GPIB”) 142; however, other suitable interfaces may be utilized. Imagesfrom the phosphorescent screen image intensifier 114 may then begathered at each phase delay by camera 112 to determine phase-sensitivefluorescence intensity. In one embodiment, each image, may contain128×128 pixels of intensity counts and be acquired by use of a 0.4-s CCDexposure time. A fast Fourier transform of the phase-sensitivefluorescence intensities acquired at each CCD pixel may be performed tocompute the amplitude and phase-shift of excitation light 108 at eachCCD pixel. Prior to generating a 3D image of fluorescent target 102based on the amplitude and phase-shift of emission light 108, thepredicted amplitude and phase-shift of emission light 108 may bedetermined by the coupled frequency-domain diffusion equations (1) and(2) discussed below; however, the amplitude and phase may be determinedby any other suitable mathematical model of light propagation. However,the solution to these coupled equations require the spatial distributionof amplitude and phase-shift of incident excitation light 104. Since itis difficult to create an excitation beam 104 of intensity-modulatedexcitation light that is spatially uniform in both amplitude andphase-shift, the spatial distribution of excitation light 104 may beexperimentally determined.

In one embodiment, system 100 may be altered to experimentally determinethe amplitude and phase-shift of incident excitation light 104 in thefollowing manner: (1) linear polarizers 144 are positioned at the outputof laser diode 116 and the input of CCD camera 112; (2) aperture of lens124 may be minimized; and (3) bandpass filter 128 and band-rejectionfilter 130 are replaced with a neutral density filter 146 to providesystem 200 illustrated in FIG. 2. Generally, re-emitted light 110includes both specularly (or singly) reflected excitation light andmultiply reflected excitation light. Since singly reflected excitationlight may be representative of incident excitation light 104, re-emittedlight 110 is compared to multiply-reflected excitation light todetermine the spatial distribution of the singly reflected excitationlight and thus the spatial distribution of incident excitation light104.

In one aspect of operation, the polarization angle of linear polarizers144 are first oriented parallel to one another to predominately detectspecularly reflected excitation light. Images of phase-sensitiveexcitation light intensity of this detect light may be determined usingthe procedure described above with regard to FIG. 1, predominantlyresulting in images representative of specularly reflected excitationlight. In order to substantially remove multiply scattered excitationlight intensities from these images, the polarization angles of linearpolarizers 144 may be adjusted to determined images of phase-sensitiveexcitation light intensity for multiply reflected light. In general,linearly polarized light incident on a tissue surface and multiplyreflected by fluorescent target 102 reemerges from the tissuedepolarized with equal probability of having a polarization angleparallel or perpendicular to incident polarization angle. Therefore, tosubstantially remove image intensities representative of multiplyreflected excitation light from the images of phase-sensitive excitationlight intensity, a second set of images of phase-sensitive excitationlight intensity representative of multiply reflected excitation lightmay be acquired after orientating the polarization angles of polarizer144 perpendicular to one another. The second set of images may besubtracted from the first set of images to provide phase-sensitiveexcitation light intensity of the specularly reflected excitation light.Afterward, a fast Fourier transform of the phase-sensitive imagesrepresentative of the specularly reflected excitation light may beperformed to determine the spatial distribution of the amplitude andphase of incident excitation light 104. Once the spatial distribution ofboth incident excitation light 104 and emission light 108 aredetermined, the mathematical expression discussed in relation to FIG. 4may be used to iteratively converge upon the optical properties offluorescent target 102.

FIG. 3 is an example flow diagram illustrating a method 300 forgenerating a 3D image from light detected on a two-dimensional sensorsurface. Method 300 is described with respect to system 100 and 200 ofFIGS. 1 and 2, respectively, but method 300 could also be used by anyother system. Moreover, system 100 and 200 may use any other suitabletechniques for performing these tasks. Thus, many of the steps in thisflow chart may take place simultaneously and/or in different orders asshown. Moreover, system 100 and 200 may use methods with additionalsteps, fewer steps, and/or different steps, so long as the methodsremain appropriate.

At a high level, method 300 illustrates two processes for generating a3D image: an emission-light process and an incident-light determinationprocess. The emission-light determination process is illustrated insteps 302 to 316, and the incident-light determination step isillustrated in steps 318 to 332. Method 300 begins at step 302 where atime-varying excitation light 104 is directed at a surface area 126 ofbody 106 to stimulate fluorescent target 102. Next, at step 304,re-emitted light 110 from body 106 is filtered by band-pass filter 128to pass a subband of re-emitted light 110 including the emissionwavelength to band-rejection filter 130 while substantially rejectingother wavelengths. Band-rejection filter 130 rejects a subband of thereceived light, including the excitation wavelength, and passes emissionlight 108 to image intensifier 114 at step 306. At step 308, imageintensifier 132 is phase-locked to laser diode 116 and amplifiesemission light 108. In order to determine phase-shift data for emissionlight 108, image intensifier 132 is stepped between 0° and 360° relativeto the phase of laser diode modulation at step 310. Next, at step 312,the amplified emission light 108 is detected at a two-dimensionalsurface area 138 of CCD camera 112 to electronically recover intensitydata of emission light 108. Images of phase-sensitive fluorescenceintensities is determined for each CCD pixel at step 314. At step 316, afast Fourier transform of the phase-sensitive fluorescence intensitiesis performed to compute spatial distribution of the fluorescenceintensity and phase-shift for each CCD pixel. In short, system 100 hasnow experimentally determined the spatial distribution of emission light108 and is ready to determine the spatial distribution of incidentexcitation light 104.

The excitation-light determination step begins at step 318 where linearpolarizers 144 are installed at the output of laser diode 116 and theinput of image intensifier 114. Next, at step 320, bandpass rejectionfilter 128 and band-rejection filter 130 are replaced with neutraldensity filter 146. In order to determine the spatial distribution forexcitation light, the polarization axis of linear polarizers 144 arealigned along the same axis at step 322. A first set of images ofphase-sensitive excitation intensities is determined at step 324, whichpredominantly represents specularly reflected excitation light. At step326, linear polarizers 144 are aligned with their polarization axissubstantially perpendicular to detect multiply scattered excitationlight. A second set of images of phase-sensitive excitation intensitiesis determined at step 328, which predominantly represents multiplyreflected excitation light. In order to substantially remove multiplyreflected excitation intensities from the first set of images, thesecond set of images is subtracted from the first set of images todetermine images of phase-sensitive excitation intensities forspecularly reflected excitation light at step 330. Since specularlyreflected excitation light is representative of incident excitationlight 104, a fast Fourier transform of the images of phase-sensitiveexcitation intensities for specularly reflected excitation lightcomputes the spatial distribution of excitation intensity andphase-shift for incident excitation light 104 at step 332. Finally, atstep 334, a 3D image of fluorescent target 102 is generated based on thespatial distribution of intensity and phase-shift for incidentexcitation light 104 and emission light 108. FIG. 4 illustrates oneembodiment of determining the 3D image corresponding to step 334 above,in accordance with the penalty/modified barrier function (PMBF) method.

FIG. 4 illustrates one embodiment of step 334 of FIG. 3 for generating a3D image of fluorescence absorption coefficients. Method 400 isdescribed with respect to system 100 and 200 of FIGS. 1 and 2,respectively, but method 400 could also be used by any other system.Moreover, system 100 and 200 may use any other suitable techniques forperforming these tasks. Thus, many of the steps in this flow chart maytake place simultaneously and/or in different orders as shown. Moreover,system 100 and 200 may use methods with additional steps, fewer steps,and/or different steps to generate a 3D image of fluorescence absorptioncoefficients.

To aid in understanding the various mathematical aspects of method 400,the following tables of selected variables is listed: {right arrow over(r)} position c speed of light S source n number of nodal points ωangular modulation frequency μ_(a) _(xi) absorption coefficient due tochromophore μ_(a) _(xf) absorption coefficient due to fluorophore μ_(a)total absorption at the emission wavelength is equal to the absorptionowing to nonfluorescent chromophores μ_(a) _(xi) and fluorescentchoromophore, μ_(a) _(mf) D optical diffusion coefficient = ⅓[μ_(a)_(x,m) + μ′_(S) _(x,m) ] μ′_(S) _(x,m) isotropic scattering coefficientμ_(a) _(x,m) total absorption coefficient Φ complex fluence φ quantumefficiency τ fluorescence lifetime Ω volume I_(AC) intensity amplitude θphase-shift K global stiffness matrices b global vector in FEM stiffnessequation containing source terms E error function N_(D) number ofdetectors N_(S) number of source g gradient vector ∇E d Newton directionα step length H Hessian matrix ∇²E l lower bound of the absorptioncoefficient owing to fluorophore u upper bound of the absorptioncoefficient owing to fluorophore L¹ Natural coordinate B(μ_(a)) activeset (lower bound) B′(μ_(a)) active set (upper bound) C(μ_(a)) freevariable set ε tolerance of lower and upper bound λ Lagrange multipliersλ^(l) Lagrange multiplies for lower bound λ^(u) Lagrange multipliers forupper bound M modified barrier penalty function h_(i(μ) _(a)) boundconstrain f(h_(i)(μ_(a))) quadratic penalty/barrier term of PMBFφ(h_(i)(μ_(a))) quadratic extrapolation method S_(i) shifting parameterη barrier parameter β parameter defines how close the extrapolation willget to singularities of logarithm term q_(i) ¹, q_(i) ², q_(i) ³coefficient of the quadratic function γ reduction factor of the barrierparameter k outer iteration number j inner iteration number γ₁ ^(k), γ₂^(k), d₁ and d₂ convergence criteria ε₁ reconstruction error

Generally, method 400 uses the modified barrier penalty function toreconstruct a 3D image of fluorescence absorption coefficients. However,prior to discussing method 400 additional information will be provided.

Diffusely propagated fluorescence is predicted by the coupled diffusionequations (1) and (2) given the optical properties of the tissue. Thenumerical solution of the diffusion equation may be needed to providesufficiently accurate estimates of the fluence on the boundary of thetissue for the image reconstruction (or inverse) problem. The coupleddiffusion equation used in optical tomography at a given modulationfrequency f(ω=2πf rad) is given by: $\begin{matrix}{{{\nabla{\cdot \left\lbrack {{D_{x}\left( \overset{\rightarrow}{r} \right)}{\nabla{\Phi_{x}\left( {\overset{\rightarrow}{r},\omega} \right)}}} \right\rbrack}} + {\left\lbrack {\frac{{\mathbb{i}}\quad\omega}{c} + {\mu_{a_{xi}}\left( \overset{\rightarrow}{r} \right)} + {\mu_{a_{xf}}\left( \overset{\rightarrow}{r} \right)}} \right\rbrack{\Phi_{x}\left( {\overset{\rightarrow}{r},\omega} \right)}}} = {S\quad{on}\quad\Omega}} & (1) \\{{{\nabla{\cdot \left\lbrack {{D_{m}\left( \overset{\rightarrow}{r} \right)}{\nabla{\Phi_{m}\left( {\overset{\rightarrow}{r},\omega} \right)}}} \right\rbrack}} + {\left\lbrack {\frac{{\mathbb{i}}\quad\omega}{c} + {\mu_{a_{m}}\left( \overset{\rightarrow}{r} \right)}} \right\rbrack{\Phi_{m}\left( {\overset{\rightarrow}{r},\omega} \right)}}} = {{\phi\mu}_{a_{xf}}\quad\frac{1}{1 - {{\mathbb{i}}\quad\omega\quad\tau}}{\Phi_{x}\left( {\overset{\rightarrow}{r},\omega} \right)}\quad{on}\quad\Omega}} & (2)\end{matrix}$Here, Φ_(x) and Φ_(m) are the AC components of the excitation andemission photon fluence (photons/cm²/sec) at position r, respectively; ωis the angular modulation frequency (rad/s); c is the speed of lightwithin the medium (cm/s); μ_(a) _(xi) is the absorption due tochromophores (cm⁻¹); and μ_(a) _(xf) is the absorption due tofluorophores (cm⁻¹). D_(x,m) is the optical diffusion coefficientequivalent to ⅓(μ_(a) _(x,m) +μ_(s) _(x,m) ′), where μ_(s) _(x,m) ′ isthe isotropic scattering coefficient (cm⁻¹) and μ_(a) _(x,m) is thetotal absorption coefficient, at the respective wavelengths. The totalabsorption at the emission wavelength, μ_(a) _(m) , is equal to theabsorption owing to nonfluorescent chromophores, μ_(a) _(xi) , andfluorescent chromophores, μ_(a) _(mf) . The right hand term of equation(2) describes the generation of fluorescence within the medium. The termφ represents the quantum efficiency of the fluorescence process, whichis defined as the probability that an excited fluorophore will decayradiatively, and τ is the fluorophore lifetime (ns). Note that thesource term couples with the solution of excitation fluence described byequation (1). The numerical solutions for the excitation and emissionfluence distributions given in equations (1) and (2), respectively, maybe obtained using the Robin boundary condition. Other suitablemathematical models of light propagation and emission generation may beused to obtain the excitation and/or emission fluence distributions.

The source term S in equation (1) depends on the measurement types. Forpoint illumination, the excitation source is represented as$\begin{matrix}{S = {\sum\limits_{i = 1}^{n_{1}}\quad{S\quad{\delta\left( {\overset{\rightarrow}{r} - {\overset{\rightarrow}{r}}_{s}} \right)}}}} & (3)\end{matrix}$where r_(s) is the positional vector of the illuminating point on thesurface of the phantom; δ is the Dirac delta function representing thesource of excitation at a single point; n₁ is the total number ofsimultaneous points of illumination, which is typically one. For areaillumination, the area of the planar source can be discretized into anumber of triangular finite elements and the source within a surfaceelement in area coordinate is represented by $\begin{matrix}{S = {\sum\limits_{i}^{3}\quad{L_{i}^{1}S_{i}}}} & (4)\end{matrix}$where the L_(i) ¹ are the natural co-ordinates for the triangle andS_(i) is the source value at each node (i). The source termcorresponding to the area illumination can be approximated as thesummation of simultaneous and closely position multiple point sources asgiven by equation (4).

The Galerkin finite element method with tetrahedral in three-dimensionsmay be used to solve the differential equations (1) and (2). Localstiffness matrices may be formed for each element and then thesematrices may be combined into two global stiffness matrices, K_(x) forthe excitation equation and K_(m) for the emission equation. Thesolutions of the differential equations are obtained by solving thecomplex equations:K _(x,m)Φ_(x,m) =b _(x,m)  (5)where b_(x,m) is the right hand side of excitation and emissionequations. The discretization process converts the continuous probleminto a problem with finite number of unknowns by expressing the unknownfield variables, which are the absorption coefficients owing tofluorophore μ_(a) _(xf) in terms linear interpolation functions withineach element. The linear functions are defined in terms of the values ofthe field variables at each node of the element. Therefore, the nodalvalues of the field variables become the unknowns to the image problem.

From the solution of Φ_(x,m), the values of the amplitude, I_(AC) _(x,m), and phase shift, θ_(x,m), can be determined from the relationship:Φ_(x,m) =I _(AC) _(x,m) exp(−iθ _(x,m))  (6)

Two finite elements meshes may be employed when solving the modifiedbarrier penalty function. The first may contain 33431 nodes with 163285tetrahedral elements (denoted below as Mesh 1) while the second, morerefined mesh may contain 38919 nodes with 193137 tetrahedral elements(defined as Mesh 2). 1205 seconds and 1738 seconds of CPU time may berequired to solve the coupled excitation and emission equations for Mesh1 and Mesh 2, respectively. All the computations in this paper may beperformed on a SUN-blade 100 UNIX workstation.

Information regarding the error function will also be provide prior todiscussing method 400. The error function considered for the imagereconstruction problems is as followed (for brevity we will use μ_(a)for μ_(a) _(xf) henceforth): $\begin{matrix}{{E\left( {\mu_{a},\omega} \right)} = {\frac{1}{2}{\sum\limits_{p = 1}^{N_{D}}\quad\left\lbrack {\left( {{\log\left( Z_{p} \right)}_{c} - {\log\left( Z_{p} \right)}_{m}} \right)_{j}\left( {{\log\left( Z_{p} \right)}_{c}^{*} - {\log\left( Z_{p} \right)}_{m}^{*}} \right)} \right\rbrack}}} & (7)\end{matrix}$where E is the error function associated with measurements, Z_(p) iscomprised of referenced fluorescent amplitude, I_(A  C_(ref_(p))),and referenced phase shift θ_(ref) _(p) measured at boundary point, p,in response to planewave illumination. Specifically the referencedmeasurement at boundary point p is given by: $\begin{matrix}{Z_{p} = {I_{A\quad C_{{ref}_{p}}}{\exp\left( {{\mathbb{i}}\quad\theta_{{ref}_{p}}} \right)}}} & (8)\end{matrix}$

The experimental values of fluorescence modulated amplitude I_(AC) _(p), may be divided by the maximum fluorescence experimental value from thecollected reflectance data set, I_(AC) _(max) in order to compute thereferenced amplitude values, I_(A  C_(ref_(p)))The experimental fluorescence phase shift, at the position where themaximum fluorescence I_(AC) _(max) was obtained, may be determined asθ_(at  I_(A  C_(max)))and subtracted from all experimental fluorescence phase shift values,θ_(p) in the collected reflectance data sets to obtained referencedphase shift values, θ_(ref) _(p) . These represent the “referenced”measurements: $\begin{matrix}{{I_{A\quad C_{{ref}_{p}}} = \frac{I_{A\quad C_{p}}}{I_{A\quad C_{\max}}}};{\theta_{{ref}_{p}} = {\theta_{p} - \theta_{{at}\quad I_{A\quad C_{\max}}}}}} & (9)\end{matrix}$The predicted fluorescence I_(A  C_(ref_(p)))  and  θ_(ref_(p))and θ_(ref) _(p) may be computed by solving the coupled diffusionequations (1 and 2) using the finite element method and may be similarlynormalized to obtain the calculated values. Other suitable methods maybe used to calculate these values. The referenced data are normalized inthis manner in order to reduce, minimize, or eliminate instrumentfunctions.

The error in equation (7) is reported for values of tissue absorptionowing to fluorescent contrast agent. The subscript c denotes the valuescalculated by the forward solution and the subscript m denotes theexperimental value. N_(D) denotes the number of detectors on thesurface. The superscript * denotes the complex conjugate of the complexfunction of the referenced measurements.

The discussion regarding the coupled diffusion equations (1) and (2) andthe error function (7) provides a basis for the modified barrier penaltyfunction, M (termed hereafter as the barrier function), which is statedas followed: $\begin{matrix}{{\min\limits_{\mu_{a}}\quad{M\left( {\mu_{a},\lambda,\eta} \right)}} = {{E\left( \mu_{a} \right)} - {\eta{\sum\limits_{i = 1}^{n}\quad{\lambda_{i}^{l}{f\left( {\mu_{a_{i}} - l} \right)}}}} + {\lambda_{i}^{u}{f\left( {u - \mu_{a_{i}}} \right)}}}} & (11)\end{matrix}$where η is a penalty/barrier term. n is the number of nodal points.λ^(l) and λ^(u) are vectors of Lagrange multipliers for the boundconstraints for the lower and the upper bounds, respectively. Fromequation (11) one can see that the simple bounds are included in thebarrier function M.

Function, f may be defined as: $\begin{matrix}{{f\left( {h_{i}\left( \mu_{a} \right)} \right)} = \left\{ \begin{matrix}{\log\left( {s_{i} + \frac{h_{i}\left( \mu_{a} \right)}{\eta}} \right)} & {{{{for}\quad i} = 1},2,\ldots\quad,{{n\quad{if}\quad{h_{i}\left( \mu_{a} \right)}} \geq {{- \beta}\quad s_{i}\eta}}} \\{\varphi\left( {h_{i}\left( \mu_{a} \right)} \right)} & {{{{for}\quad i} = 1},2,\ldots\quad,{{n\quad{if}\quad{h_{i}\left( \mu_{a} \right)}} < {{- \beta}\quad s_{i}\eta}}}\end{matrix} \right.} & (12) \\{where} & \quad \\{{\varphi\left( {h_{i}\left( \mu_{a} \right)} \right)} = {{\frac{1}{2}{q_{i}^{2}\left( {h_{i}\left( \mu_{a} \right)} \right)}^{2}} + {q_{i}^{2}\left( {h_{i}\left( \mu_{a} \right)} \right)} + q_{i}^{3}}} & (13)\end{matrix}$h_(i)(μ_(a)) represents all the bound constraints {(μ_(a) _(i) −l) and(u−μ_(a) _(i) )}; and s_(i) is positive-valued shifts used in theenlargement of the feasibility region for the lower and the upperbounds, respectively. In one embodiment, the s_(i) is set to 1.0 for allcomputations. The coefficients q_(i) ¹, q_(i) ², and q_(i) ³ of thequadratic function in equation (13) are chosen so that their values, andtheir first and second derivatives match those of the logarithmic termin equation (12) at the point h_(i)(μ_(a))=−βs_(i)η. The coefficients ofthe quadratic as given as: $\begin{matrix}{q_{i}^{1} = \frac{- 1}{\left( {s_{i}{\eta\left( {1 - \beta} \right)}} \right)^{2}}} & (14) \\{q_{i}^{2} = \frac{1 - {2\beta}}{s_{i}{\eta\left( {1 - \beta} \right)}^{2}}} & (15) \\{q_{i}^{3} = {\frac{\beta\left( {2 - {3\beta}} \right)}{\left( {2\left( {1 - \beta} \right)^{2}} \right)} + {\ln\left( {s_{i}\left( {1 - \beta} \right)} \right)}}} & (16)\end{matrix}$The coefficients are dependent on both β and η, and thus have to beupdated for any changes of β and η. The values of β defines how close tothe singularities the logarithmic terms described in equation (12) areextrapolated. Consequently we define a very simple differentiablemodified barrier penalty function by equation (11) that is used toassess the quality of the solution. This function has an extremelysimple structure so that the overhead for using the barrier functioninstead of the original error function given by equation (7) isnegligible for function evaluation or for computing gradients andHessian as discussed below.

The PMBF method 400 is performed in two stages in, an inner and an outerstage illustrated in FIG. 4. The outer iteration updates the Lagrangemultiplier λ, parameter β, and the barrier parameter η. Formulationsused to update these parameters at each outer iteration are provided inthe discussion below. Using the calculated values of the parameter λ, anunconstrained technique may be applied to minimize the penalty barrierfunction described by equation (11). Once satisfactory convergencewithin the inner iteration is reached, the variables in the outeriteration are recalculated to check for satisfactory convergence in theouter iteration. If unsatisfactory, the Lagrange multipliers, thebarrier parameters are updated and another unconstrained optimization isstarted within the inner iteration.

Turning to method 400, method 400 begins at step 402 where initialvalues of μ_(a), η and β are selected. The initial value for μ_(a), theabsorption coefficient due to fluorophore, is a positive value typicallyselected from the range of 0 to 0.1 cm⁻¹. The initial value for η, thebarrier parameter, is an arbitrary positive value and typically 10⁻¹ or10⁻². The initial value for β is a parameter that defines how close tothe singularities the logarithmic terms in the modified barrier penaltyfunction are extrapolated. The values of β are typically not too small(e.g., 0≦β≦0.5) or very close to one (i.e., β=0.99). For example, thevalue of β may be set to 0.9. Next, at decisional step 404, if theinitial values of λ^(l), an initial Lagrange multiplier for a lowerbound, and λ^(u), an initial Lagrange multiplier for an upper bound, areselected, then, at step 406, the values are set to equal values such as,for example, 1, 10, 100 and 1000. If the initial values of λ^(l) andλ^(u) are not initially selected, then, at step 408, the values of λ^(l)and λ^(u) are determined. In one embodiment, the initial values of λ^(l)and λ^(u) are determined using a vector λ⁰ that minimizes theleast-squares unconstrained problem: $\begin{matrix}{\min\limits_{\lambda \geq 0}{{{\nabla_{\mu_{a}}{E\left( \mu_{a}^{0} \right)}} - {\lambda\quad{\nabla_{\mu_{a}}{f\left( {h\left( \mu_{a}^{0} \right)} \right)}}}}}_{2}^{2}} & (17)\end{matrix}$this gives the formula:λ(μ_(a) ⁰)=A(μ_(a) ⁰)∇_(μ) _(a) E(μ_(a) ⁰)∇_(μ) _(a) f(h(μ_(a) ⁰))  (18)whereA(μ_(a) ⁰)=(∇_(μ) _(a) f ^(T)(h(μ_(a) ⁰))∇_(μ) _(a) f(h(μ_(a)⁰)))⁻¹  (19)The formulation of equation (18) requires calculation of all constraintgradients and the gradients of the error function E(μ_(a) ⁰) andinversion of the matrix A(μ_(a) ⁰). The inversion of this large (n×n)matrix in the reconstruction problem can be very costly. To reduce thecomputational cost, only the diagonal terms of the equation (18) may beused. Hence there is no need to invert the matrix. Using thisleast-square procedure the initial Lagrange multipliers for imagereconstruction problem may be computed. After the initial values ofμ_(a), η, β, λ^(l), and λ^(u) are set, the indices k, which identifiesthe current step of the outer stage, is set to zero at step 410.

Method 400 now executes the inner stage of the PMBF method by settingthe indices j, which represents the current iteration of the innerstage, to zero at step 412. Next, at step 414, the barrier function, thegradient of the barrier function and the Hessian of the barrier functionare calculated. The barrier function may be calculated by equation (11).During this calculation, an unconstrained optimization method may beused to optimize the unconstrained barrier function of equation (11) foreach set of parameters λ^(k) and η^(k). The truncated Newton method withtrust region (TN) may be used for the unconstrained optimizationproblem. The gradient based truncated Newton method is efficient andsuitable for large-scale problem. The gradient of the barrier functionmay given by $\begin{matrix}{{\nabla_{\mu_{a}}{M\left( {\mu_{a},\lambda^{k},\eta^{k}} \right)}} = \left\{ \begin{matrix}{{{\nabla_{\mu_{a}}{E\left( \mu_{a} \right)}} - {\eta^{k}{\sum\limits_{i = 1}^{n}\quad{\frac{\lambda^{k}}{\left( {\frac{h_{i}\left( \mu_{a} \right)}{\eta^{k}} + s_{i}} \right)}{\nabla{h_{ix}\left( \mu_{a} \right)}}\quad{if}\quad{h_{i}\left( \mu_{a} \right)}}}}} \geq {{- \beta^{k}}s_{i}\eta^{k}}} \\{{{{\nabla_{\mu_{a}}{E\left( \mu_{a} \right)}} - {\eta^{k}{\sum\limits_{i = 1}^{n}{\lambda^{k}\left\{ {{q_{i}^{1}\left( {h_{i}\left( \mu_{a} \right)} \right)} + q_{i}^{2}} \right\}\quad{if}\quad{h_{i}\left( \mu_{a} \right)}}}}} < {{- \beta^{k}}s_{i}\eta^{k}}}\quad}\end{matrix} \right.} & (20)\end{matrix}$The gradients of the error function ∇_(μ) _(a) E may be calculated bythe analytical method using reverse differentiation technique. TheHessian of the barrier function may be determined as followed:$\begin{matrix}{{\nabla_{\mu_{a}}^{2}{M\left( {\mu_{a},\lambda^{k},\eta^{k}} \right)}} = \left\{ \begin{matrix}{{\nabla_{\mu_{a}}^{2}{E\left( \mu_{a} \right)}} - {\eta^{k}{\sum\limits_{i = 1}^{n}\frac{\lambda^{k}}{\left( {{h_{i}\left( \mu_{a} \right)} + {s_{i}\eta^{k}}} \right)^{2}}}}} & {{h_{i}\left( \mu_{a} \right)} \geq {{- \beta^{k}}s_{i}\eta^{k}}} \\{{\nabla_{\mu_{a}}^{2}{E\left( \mu_{a} \right)}} - {\eta^{k}{\sum\limits_{i = 1}^{n}{\lambda^{k}\left\{ q_{i}^{1} \right\}}}}} & {{h_{1}\left( \mu_{a} \right)} < {{- \beta^{k}}s_{i}\eta^{k}}}\end{matrix} \right.} & (21)\end{matrix}$∇_(μ) _(a) ²E(μ_(a)).d maybe calculated by the finite difference scheme,$\begin{matrix}{{{H\left( \mu_{a} \right)} \cdot d} = {\frac{1}{\delta}\left\{ {\left\lbrack {g\left( {\mu_{a} + {\delta\quad d}} \right)} \right\rbrack - {g\left( \mu_{a} \right)}} \right\}}} & (22)\end{matrix}$where ∇_(μ) _(a) E(μ_(a))=g(μ_(a)) is the gradient; ∇_(μ) _(a)²E(μ_(a))=H(μ_(a)) is the Hessian of the error function E(μ_(a)); d isthe Newton direction; and δ is (approximately) the square root of themachine precision. It may not safe to perform finite differencingdirectly to ∇_(μ) _(a) ²M because of the singularity of the logarithmicfunction. The final summation in the formulation for ∇_(μ) _(a) ²M maybe computed straightforwardly from the formulas above.

If the absolute value of the gradient of the barrier function is lessthan 10⁻⁵ at decisional step 416, then the inner stage ends andexecution proceeds to step 430 of the outer stage. Otherwise, ∇_(μ) _(a)²M(μ_(a) ^(j))d^(j)=∇_(μ) _(a) M(μ_(a) ^(j)) is solved at step 418 todetermine the Newton direction d. The Newton direction d is computed asan approximation solution of the Newton equations:(∇² M)d=−∇M  (23)where ∇²M≡∇²M(μ_(a)) is the Hessian matrix of second derivatives at thecurrent solution designated by a vector μ_(a). The approximate solutionof equation (23) may be obtained by the linear conjugate-gradient methodbut the solution may alternatively be determined by any suitable method.This iterative method is “truncated” before the exact solution isobtained. The conjugate-gradient method (CG) corresponds to minimizingthe quadratic model $\begin{matrix}{{\min\limits_{d}{Q(d)}} = {{\frac{1}{2}d^{T}{\nabla^{2}{Md}}} + {d^{T}{\nabla\quad M}}}} & (24) \\{{{s.t.\quad l} - \mu_{a}} \leq d \leq {u - \mu_{a}}} & (25)\end{matrix}$as a function of d over a sequence of subspaces of increasing dimension.Criteria must be determined and applied to find out whether convergencehas been achieved. This method has two levels of iterations, with the‘outer-level’ iteration determining the step length, and the‘inner-level’ iteration computing the Newton direction using conjugategradient method. A further issue for unconstrained minimization is whatto do when the current Hessian ∇²M is not positive definite. Severalapproaches are used, among the most popular of which are so-calledmodified Newton method in which d_(j) satisfies the linear system,(∇² M(μ_(a) ^(j))+R)d ^(j) =−∇M  (26)where R is a positive semidefinite matrix. Different techniques haveused to calculate R so that (∇²M(x^(k))+R) is positive definite. Analternative globalization strategy in unconstrained optimizationembodied in trust region methods, is to choose the step d^(j) tominimize the local Taylor-series model subject to a limitation on thesize of ∥d∥, where various norms may appear in the trust-regionconstraints. To alleviate this difficulty, the conjugate gradient methodwith trust region may be used for unconstrained optimization problems.The procedure is described in the following paragraph as followed:

CG methods generate sequences p^(i), y^(i), where p^(i) approximatesiteratively the Newton direction d^(j) and the conjugate directiony^(i). The update scheme isp ^(i+1) =p ^(i)+Λ^(i) y ^(i)  (27)y ^(i+1) =r ^(i+1)+ρ^(i) y ^(i)  (28)where Λ^(i) minimizes the quadratic function Q(d) along the directiony^(i), and ρ^(i) is chosen so as to maintain conjugacy between thedirections y^(i). The CG terminates when the residual r^(i)=−∇Q(p^(i))is in norm below a prescribed tolerance or when a negative curvaturedirection τ₁=(s^(i))^(T)∇²Ms^(i)≦0 is encountered. This algorithm doesnot stop, but computes τ₁>0 so that ∥p^(i)+τ₁y^(i)∥=Δ, where Δ is thetrust region. Then update rule is d^(i+1)=d^(i)+τ₁y^(i). The algorithmstops if the residual ∥r^(i)∥ is below a given tolerance, or if(y^(i))^(T)∇²My^(i)=0.

Next, at step 420, a new estimate of the absorption owing to fluorophoreis determined by μ_(a) ^(j+1)=μ_(a) ^(j)+αd^(j), where α is the steplength. The step length may be determined by an Armijo line searchalgorithm as disclosed in M. S. Bazaraa, H. D. Sherali, and C. M.Shetty, Nonlinear programming theory and algorithms, John Wiley & SonsInc, New York., 1993. The Newton direction d may be required to satisfyd^(T)∇M<0 (i.e., it is a descent direction for M at the point μ_(a)).The step length α>0 may be chosen to guarantee that M({circumflex over(μ)}_(a))<M(μ_(a)), along with Wolfe's conditions disclosed in P. Wolfe,“Convergence condition for ascent method,” SIAM Rev., 11, 226-253, 1969,which are designed to guarantee convergence to a local minimum.

At decisional step 422, its is determined whether μ_(a) ^(j+1) is anactive variable. In general, if the variables of an unconstrainedproblem are very close to a bound in the inner iteration, then thesevariables are fixed (call active variables). Hence the start of the nextunconstrained optimization uses the free variables only. The variablesin B(μ_(a)) and B′(μ_(a)) are called active variables and the variablesin C(μ_(a)) are called free variablesB(μ_(a))={i:l≦μ _(a) _(i) ≦l+ε}  (25)B′(μ_(a))={i:u−ε≦μ _(a) _(i) ≦u}  (26)C(μ_(a))={i:l+ε≦μ _(a) _(i) ≦u−ε}  (27)where i is the number of nodal points i=1, . . . n, l is the lower boundof the absorption coefficient owing to fluorophore, u is the upper boundof the absorption coefficient owing to fluorophore, and the tolerance Eshould be sufficiently small, such thatl+ε<u−ε  (28)For example, the tolerance ε may be set to 10⁻³, 10⁻⁴, 10⁻⁶, or othersuitable values. All variables in the active set are fixed during theiteration, while the other variables in set C are free. After eachiteration, the active variables may be updated from current freevariables, i.e., variables may be found that lie within the bounds andadd to previous active variables. The active sets B and B′ are updatedat the inner iteration and the process is continued until the gradientof the barrier function M described by (11) can not be reduced furtherwith the current free variables.

The variables μ_(a) _(i) in C={μ_(a) _(i) |l≦μ_(a) _(i) ≦u} have threemutually exclusive possibilities: $\begin{matrix}{{(i)\quad 0} \leq {\mu_{a_{i}} - l} \leq {\nabla_{\mu_{a_{i}}}{M\left( {\mu_{a},\lambda,\eta} \right)}}} & (29) \\{{({ii})\quad{\nabla_{\mu_{a_{i}}}{M\left( {\mu_{a},\lambda,\eta} \right)}}} \leq {\mu_{a_{i}} - u} \leq 0} & (30) \\{{({iii})\quad\mu_{{a_{i} - u} \leq}{\nabla_{\mu_{a_{i}}}{M\left( {\mu_{a},\lambda,\eta} \right)}}} \leq {\mu_{a_{i}} - l}} & (31)\end{matrix}$The vector of μ_(a) _(i) that satisfies (i) or (ii) are now thedominated variables. The ones satisfying (i) are said to be dominatedabove; and those satisfying (ii) dominated below; and those satisfying(iii) are the floating variables.

If a convergence condition is satisfied at decisional step 424, theninner stage ends and method 400 proceeds to step 430 of the outer stage.In the illustrated embodiment, the inner iteration can be consideredfully converged if either of the following two conditions are satisfied.The first convergence condition is $\begin{matrix}{{\frac{\nabla_{\mu_{a}}{M\left( {\mu_{a}^{j},\lambda^{k},\eta^{k}} \right)}}{1.0 + {\mu_{a}^{j}}}} < \tau_{1}^{k}} & (32)\end{matrix}$where k is the number of the subproblem and j is the number of inneriteration. If τ₁ ^(k) is less than 10⁻⁵ then the inner iteration isstopped and the outer iteration are continued. The second convergencecondition is satisfied if the fractional change d₂ given by (33) of thebarrier function M of equation (11) of consecutive inner iterations isless than ε₀=10⁻⁵ then the inner iteration is stopped and the outeriteration are continued. $\begin{matrix}{d_{2} = \frac{{{M\left( {\mu_{a}^{j + 1},\lambda^{k},\eta^{k}} \right)} - {M\left( {\mu_{a}^{j},\lambda^{k},\eta^{k}} \right)}}}{1.0 + {{M\left( {\mu_{a}^{j},\lambda^{k},\eta^{k}} \right)}}}} & (33)\end{matrix}$If neither convergence condition is satisfied at decisional step 424,then, at step 426, the inner iteration indices j is set to j+1. If theinner iteration indices j does not satisfies the specified number ofinner iterations (e.g., 5, 10, etc) at decisional step 428, the inneriteration returns to step 418. Otherwise, the inner iteration ends, andexecution returns to the outer iteration at step 430.

At step 430, the outer iteration indices k is set to k+1. If the outeriteration indices k satisfies the specified number of outer iterations(e.g., 6) at decisional step 432, then, at step 434, the absorptioncoefficient due to fluorophores μ_(a) is set to the current estimatedvalue μ_(a) ^(j+1) and execution ends. If the outer iteration indices kdoes not satisfy the specified number of outer iterations at decisionalstep 432, then execution proceeds to decision step 436. If a convergencecriteria is satisfied at decisional step 436, then execution ends. Theouter iteration of the PMBF can be considered fully converged if itsatisfies any of the following conditions: $\begin{matrix}{\gamma_{1}^{k + 1} = {\frac{{{\nabla_{\mu_{a}}{E\left( \mu_{a}^{k + 1} \right)}} - {\eta^{k}{\sum\limits_{i = 1}^{n}{\lambda_{i}^{k}{\nabla_{{\mu\quad}_{a}}{h_{i}\left( \mu_{a}^{k + 1} \right)}}}}}}}{1.0 + {\mu_{a}^{k + 1}}} < \varsigma}} & (34) \\{\gamma_{2}^{k + 1} = {\frac{{\eta^{k}{\sum\limits_{i = 1}^{n}{\lambda_{i}^{k}{\nabla_{\mu_{a}}{h_{i}\left( \mu_{a}^{k + 1} \right)}}}}}}{1.0 + {\mu_{a}^{k + 1}}} < \varsigma}} & (35)\end{matrix}$where k is the number of subproblem. If ç is less than 10⁻⁵, then theinner iteration is stopped. If the fractional change d₁ given byequation (36) of the barrier function M given by equation (11) for theconsecutive outer iterations is less than ε₀=10⁻⁵, then the outeriteration is stopped. $\begin{matrix}{d_{1} = \frac{{{M\left( {\mu_{a}^{k + 1},\lambda^{k + 1},\eta^{k + 1}} \right)} - {M\left( {\mu_{a}^{k},\lambda^{k},\eta^{k}} \right)}}}{1.0 + {{M\left( {\mu_{a}^{k},\lambda^{k},\eta^{k}} \right)}}}} & (36)\end{matrix}$If the outer iteration convergence criteria are not satisfied atdecisional step 436, then, at step 438, λ^(l), λ^(u), and β are updated.The Lagrangian multipliers may be updated as followed: $\begin{matrix}{\lambda^{k + 1} = \left\{ \begin{matrix}{\eta^{k}{\lambda_{i}^{k}\left( \frac{1}{{\eta^{k}s_{i}} + {h_{i}\left( \mu_{a} \right)}} \right)}} & {{{{for}\quad i} = 1},2,\ldots\quad,{{n\quad{h_{i}\left( \mu_{a} \right)}} \geq {{- \beta^{k}}s_{i}\eta^{k}}}} \\{\eta^{k}{\lambda_{i}^{k}\left( {{q_{i}^{1}{h_{i}\left( \mu_{a} \right)}} + q_{i}^{2}} \right)}} & {{{{for}\quad i} = 1},2,\ldots\quad,{{n\quad{h_{i}\left( \mu_{a} \right)}} < {{- \beta^{k}}s_{i}\eta^{k}}}}\end{matrix} \right.} & (37)\end{matrix}$where λ_(i) represents all the Lagrange multipliers for the boundconstraints (λ^(l) and λ^(u)) and h_(i) (μ_(a)) represents all the boundconstraints. η may be updated by a simple reduction rule:η^(k+1)=η^(k)/γ  (38)where γ>1 is defined by the user, and normally γ=2 or γ=10, and k is thenumber of the outer iteration.

Next, at step 440, μ_(a) _(i) is updated at nodal point i, after k^(th)outer iteration as follows: $\begin{matrix}{\left( {\hat{\mu}}_{a}^{k} \right) = \left\{ \begin{matrix}l & {{{if}\quad 0} \leq {\left( \mu_{a}^{k - 1} \right)_{i} - l} \leq \left( {\nabla_{\mu_{a}}M^{k - 1}} \right)_{i}} \\u & {{{if}\quad\left( {\nabla_{\mu_{a}}M^{k - 1}} \right)_{i}} \leq {\left( \mu_{a}^{k - 1} \right)_{i} - l} \leq 0} \\\left( \mu_{a}^{k - 1} \right)_{i} & {otherwise}\end{matrix} \right.} & (39)\end{matrix}$Execution returns to step 412 to return to the inner stage.

FIGS. 5A-C illustrate x-z views of a 3D distribution of a fluorescenttarget which provides example data generated according to the teachingsof the invention. These reconstructions used the following initialvalues μ_(a) _(xf) ⁰=0.03 cm⁻¹ and λ⁰=1000, μ_(a) _(xf) ⁰=0.003 cm⁻¹ andλ⁰=1000, and μ_(a) _(xf) ⁰=0.003 cm⁻¹ and λ⁰=1000, respectively. Atissue phantom consisting of an 8×8×8 cm³ cube filled with 1% lipidsolution to mimic the scattering properties of tissue was used in thesystem of FIGS. 1 and 2. The target to be reconstructed in these initialstudies was a 1 cm³ cube filled with 1% lipid and 1 μM Indocyanine Greensolution placed 1 cm from the imaging surface. The top surface of thetissue phantom was illuminated with an expanded beam of 785-nmexcitation light which was intensity modulated at 100 MHz. The surfacearea illuminated was approximately 9 cm². The emitted fluorescent signalin this area was imaged using a gain modulated intensified chargecoupled device (ICCD) operated under homodyne conditions to enabledetermination of the amplitude and phase-delay of the light emitted onthe top of the surface. Using a combination crossed and parallelpolarizers on the incident and collected light, the spatial distributionof the amplitude and phase delay of the incident excitation light wasdetermined. Next upon placing an 830 nm interference and imagingholographic filters at the photocathode of the intensifier, theamplitude and phase delay of the re-emitted fluorescent light owing tothe submerged target was determined across the surface of the phantom.

One advantage of this procedure is that optical fibers are not requiredto illuminate or collect the light signal unlike the traditional pointillumination/collection measurement schemes used in optical tomography.An additional advantage of this procedure is the ability to recontruct3D images from data collected at a 2D sensor surface.

Although the present invention has been described in detail, variouschanges and modifications may be suggested to one skilled in the art. Itis intended that the present invention encompass such changes andmodifications as falling within the scope of the appended claims.

1. A method of imaging, comprising: directing time-varying excitationlight at a surface area of a tissue, the tissue comprising a fluorescenttarget; modulating the intensity of the excitation light at a radiofrequency; detecting, substantially at a two-dimensional sensor surface,time-varying emission light from the fluorescent target in response tothe time-varying excitation light stimulating the fluorescent target;filtering the time-varying emission light to reject excitation lightre-emitted from the material; determining a spatial distribution ofamplitude and phase-delay of the detected time-varying emission light;detecting time-varying excitation light specularly reflected by thefluorescent target; determining a spatial distribution of amplitude andphase-delay of the time-varying excitation light incident on the surfacearea based on the reflected time-varying excitation light; andgenerating a three-dimensional image of the fluorescent target based onthe spatial distributions of the incident time-varying excitation lightand the time-varying emission light.
 2. A method of imaging, comprising:directing time-varying excitation light at a surface area of a lightscattering material, the material comprising a fluorescent target;detecting, substantially at a two-dimensional sensor surface,time-varying emission light from the fluorescent target in response tothe time-varying excitation light stimulating the fluorescent target;filtering the time-varying emission light to reject excitation lightre-emitted from the material; and generating a three-dimensional imageof the fluorescent target based on the detection substantially at thesensor surface.
 3. The method of claim 2, wherein the fluorescent agentcomprises indocyanine green (“ICG”).
 4. The method of claim 2, whereinthe surface area is greater than 1 millimeter squared (mm²).
 5. Themethod of claim 2, further comprising determining a spatial distributionof amplitude and phase-delay of the detected time-varying emissionlight.
 6. The method of claim 5, further comprising: detectingtime-varying excitation light specularly reflected by the fluorescenttarget; and determining a spatial distribution of amplitude andphase-delay of the time-varying excitation light incident on the surfacearea based on the reflected time-varying excitation light.
 7. The methodof claim 6, wherein in the three dimensional image of the fluorescenttarget is based on the spatial distributions of the incidenttime-varying excitation light and the time-varying emission light. 8.The method of claim 2, wherein the time-varying excitation lightcomprises a planar wave.
 9. The method of claim 2, wherein the materialcomprises tissue.
 10. The method of claim 2, further comprisingmodulating the intensity of the excitation light at a radio frequency.11. The method of claim 6, wherein detecting time-varying excitationlight singularly reflected by the fluorescent target comprises:detecting excitation light emitted from the light scattering material;detecting multiply scattered excitation light emitted from the lightscattering material; and subtracting the detected multiply scatteredexcitation light from the detected excitation light emitted from thelight scattering material to determine the specularly reflectedexcitation light.
 12. The method of claim 2, wherein the two-dimensionalsensor surface is approximately 16 cm².
 13. A system for imaging,comprising: a laser diode operable to direct time-varying excitationlight at a surface area of a light scattering material, the materialcomprising a fluorescent target; an image intensifier operable todetecting, substantially at a two-dimensional sensor surface,time-varying emission light from the fluorescent target in response tothe time-varying excitation light stimulating the fluorescent target;one or more optical filters operable to filter the time-varying emissionlight to reject excitation light re-emitted from the material; and animaging apparatus operable to generate a three-dimensional image of thefluorescent target based on the detection substantially at the sensorsurface.
 14. The system of claim 13, wherein the fluorescent agentcomprises indocyanine green (“ICG”).
 15. The system of claim 13, whereinthe surface area is approximately 9 centimeters squared (cm²).
 16. Thesystem of claim 13, the imaging apparatus further operable to determinea spatial distribution of amplitude and phase-delay of the detectedtime-varying emission light.
 17. The system of claim 13, furthercomprising a plurality of linear polarizers operable to pass reflectedexcitation light in a first configuration and multiply-reflectedexcitation light in a second configuration; and the image apparatusfurther operable to determine a spatial distribution of amplitude andphase-delay of the time-varying excitation light incident on the surfacearea based on the reflected and multiply reflected excitation light. 18.The system of claim 17, wherein in the three dimensional image of thefluorescent target is based on the spatial distributions of the incidenttime-varying excitation light and the time-varying emission light. 19.The system of claim 13, wherein the time-varying excitation lightcomprises a planar wave.
 20. The system of claim 13, wherein thematerial comprises tissue.
 21. The system of claim 13, wherein anintensity of the time-varying excitation light is modulated at a radiofrequency.